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PREFACE 


This manual presents a brief discussion of the theory and 
a detailed description of the computer program for an airfoil in- 
version technique developed by Eppler. The program represents re- 
vision of a deck supplied to us by Dr. Stanley Miley. Contribu- 
tions by Dr. T. E. Edwards and by Mr. R. H. Awker are also grate- 
fully acknowledged. While the program has been set up specifi- 
cally for the University of Illinois IBM 360/75 computer, its 
modification to other computers should be straightforward. 

This work was funded by NASA under grant NGR 14-005-144 as 
part of a low speed airfoil study. 



I. THEORY 


1 2 

The method developed by Eppler ^ is an inverse conformal mapping technique 
that determines the x and y coordinates from a given velocity distribution. The 
two planes involved are shown in figure 1. The C plane shows the flow about a 
circular cylinder, while the z plane represents the flow about the airfoil. The 
velocity in the z plane is given in terms of coordinates determined in the ? 
plane, z and ? are defined as: 


z - X * iy 

CO 

? = ^ -f iri = re^^ 

C2) 


The flow in the S plane is such that the rear stagnation point falls on the 
real axis at C = 1. 

There exists a transformation of the C plane to the z plane such that the 
z plane represents parallel flow about a closed airfoil at an angle of attack a. 
Since C 1 represents a stagnation point, the Kutta condition requires that this 
must transform to the trailing edge of the airfoil. As this is to be an infinite 
parallel flow, z (“) = ” ^^st be real. The general function that satis- 

fies these requirements is 


z 


CO 


CO = S 

v=o 


3 



where Sj is real, but not equal to zero. 

The complex potential in the ? plane can be represented as 

F CO = $ + i'i' = + ^3 - ^ In ? 

where F is given by 

F Si 4irC sina 


The complex velocity in the z plane is given by 


V -l.dPl-l dF/del 

Uzl |dz/dC| 




C33 


C4} 

CS) 

C63 


-~yr 
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The inverse o£ this will be used, or 


dz _ dz/dC 
dF " dF/dS 


( 7 } 


In order to prevent an undesirable root, this can be written as 


, dz , dz Yt, 

dF dz ~ ^ ds 


( 8 ) 


The velocity vector in the z plane can be introduced as V = 


Ve^®. Then 


§- Ve-« 


C 9 ) 


Therefore, 


In ~ ” In V + i0. 


CIO) 


The real part o£ In || is then -In V. Outside the boundary of the unit circle. 


dz . 


In is regular, and can be calculated when the real part is known on the 
boundary. Since F C?) is also known (eq (4) 3 , equation C8) can be solved for 
. z (z) then can be found by integration, z (C) must be of the form of 
equation C3) to match the boundary conditions. Therefore, the problem is to 
find an equation for In which results in a satisfactory . 

From equations C4) and (B) 


dF ia „ , 1 


- - - 1 ^ “2ia. 

1 ) C + e ) • 


( 11 ) 


In light of the singularities involved at the stagnation points. In 


dz 

dF 


can be represented by 


In 


dz 

dF 


ia ,1 


-2ia. 


= - In G - In C~ + e + E ib„) ? 


.-m 


m-o 


n' 


( 12 ) 


Using this equation and equations (8) and (11) , 


in ^ = in (1 - r Ca„ ^ ib J c'” 

^ ^ m=o 


( 13 ) 
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This results in 


S (a_‘+ ib ) c 

II15S0 


Expanding tbis yields 


dz ,, 1, >o * 

gj- = Cl - jD e e 

t 

z (?) must be of tbe form of equation (3) > so ^ must be of the form 


if = 


If we let '*■ equation (15) can be further expanded to 


dz 1 A 2 A ® *1 *1^ ^1^ ^2 


2 1 3 J 


+•••)■ 


Comparing equations (16) and (17) yields 

A ^ A ^ A 

6i = i*^*2T*3T*-" = 


o=z 

Since is real, is also real. At infinity, we want ^ indicating that 

the flow at infinity in the circle plane is in the same direction as the flow at 

dz 

infinity in the airfoil plane. From equation (16) , ^ at infinity. There- 


fore, 


3j - ."c = 1 


Because of this, A =0. 

o 

1 

Congsaring the — term yields 


Pi A, 3^ • 
i. + 1 = 0 

? — — ^ 


: ;-5-' 

So a, = 1, = 0. 

At infinity, we can arbitrarily set the velocity equal to unity. From 
equation (12) ^ in the limit as dz" ” G=1 and In C=0. 

The problem that remains is to define the a^^^ and b^^^ that have not yet 
been defined such that, along the surface of the airfoil, the Velocity assumes 
the prescrib ed values . 

Along the surface of the circle plane, C = Using this, equation (12) 

results in 

In ^ = -In £e^“ (e^'*’ + e^^“)] ♦ Z C\ + %) e C21) 

m=o 

If we use the substitutions 

C^3 = Ca^ cos m^ + bj^ sin m^) (22) 

m=o 

CO 

Q (d)} = S C- a sin m6 + b„ cos mi) (23) 

lit IE 

m=o 

1 , 

equation (21) can be written as - 

In § = ln[n^'*'''2 * P C« 

+ iQ(<{)) (24) 

This can be written as 

-In V((|)) + i0((l)) = - In 2 1 cos C<{)/2-a) |+P (<|)) + i[|- + Q((^3) 

, - {{v}}] (25) 

% where {{ir}} is given as 



U-:tU - ^ C0<(f)<v + 2a) 
tiTTj-j - ^ ('n-+2a<(j)<2Tr) 


( 26 ) 


The {{it}} texm is necessaxy due to the shift in the direction of the velocity 
at the stagnation point. 

The real part of equation (25) can be rearranged into the form 


P C<ij) = In 2jcos c| - cO I -In V ((^). 


(27) 


Through harmonic analysis, the nnd b^*s can be determined from 
equation (27). However, we must have a^ = 0, - 1, and = 0, due to 

equations (19) and (20). Therefore, 


/'^P(<|>) = 0 


(28) 


2it. 


J P((|)) COS(J» d^ =: TT 


(29) 


2tt. 


/"'"P(^) sintj. d<|) = 0 


(30) 


The velocity distribution we specify must meet these requirements. 
By integrating equation (14) , the transformation 


S Ca„ * ib J r” 

ZK) = /Cl- 1/53 e , d? 


(31) 


M A, 


can be derived. This yields the flow for the entire z plane. If 5 » e ^ is 
entered into this transformation, the resulting z = x + iy will yield the profil 
of the airfoil. The results are 


x(4)) - /-4 sin ^ 1 

cos (^ - 

a). 1 

1 V(cji) 

c| t QC*33 d* 

(32) 

y((f)) - /-4 sin “ 1 

r4> 

COS Cy - 

a) 

1 v(^) 

c|-+ ClCWD d'l>. 

(33) 


The only quantity remaining to be defined, then, is QC4')- However Q(4>) 
is a conjugate harmonic function of PCiji), and can be derived from the formula 
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2tt 

QW) / PCa) cot Cr^3 da 

o 

Given a velocity distribution that yields a P(ij)} such that equations C28) 
through (30) are satisfied, and an angle of attack, the x and y coordinates of 
the desired airfoil can be generated using equations (32) and (33). The angle 
of attack, a, need not be held constant, but can be a function of Thus, 
the upper surface can be designed at a different (higher for reasonable airfoils) 
angle of attack than the lower surface, or even different portions of the upper 
(or lower) surface can be designed at different angles of attack. 

The profile of the airfoil is determined by a and b„. Therefore, for 
^ ■'mm 

a fixed profile, and are fixed. Altering the angle of attack will not 
alter the airfoil profile, and, therefore, will not alter or bj^^. This means 
that P(4>) is independent of the angle of attack. Equation (27) can, then, be 
written in the form of 


P(4>) == In 2jcos ( f ct (<Ji))| " In V* (<fi) 


(3S) 


where V (ij>) is the specified velocity at the point on the airfoil corresponding 

•fc 

to and a (^) is the corresponding angle of attack. At any angle of attack, a, 
then, the velocity can be given as 


V(<|>,a) 


cos (4)/ 2 -g) 

cos ((j)/2-a* (4 j)) 


V* ((j)) 


(36) 


The circle plane can be divided into I segments, as in figure (2), where 

a 

indicates the stagnation point. The angle of 

1» ■“’a L 

attack specification takes the form 

a (4') = = constant for (37) 

and the velocity takes the form 


V m = V^W^(c^) 


C38) 



_Su 


where is a constant for ^ given as 

w»3 = [l+K '”fl - 0.36{ {‘=f* 

on the upper surface. On the lower surface, the velocity distribution is similar, 
but different values of Kpj, y, and 4»g S'Ve used, indicated by Kj^, y, and 
respectively. The terms in the double brackets of equation .(38) are defined by 

{f((/i) fCt|))> 0} 

{{f(«»={ } (40) 

{ 0 fC<J))< 0} 

Equation (39) can be considered to be of the form 
W(4>) - W^(({,) Wg(^) 
or, on the lower surface 

w.CW'VW'ifsCW " C42) 

The W (or W ) term produces the major pressure recovery. Tliis is in the 
w w 

form specified by Wortmann,^ in which the shape parameter is held constant, which 
delays separation. The W (or W ) terra develops the cusp distribution. It is 
generally applied to the last 3-5% of the airfoil length. Outside of the range 
of the specified region, ^><jj for W or ^>(|> for W (or ^<i|) or (Jj on the lower 

V/ W <5 S. n S 

surface), the value of W ^ or W is one, as dictated by equation (40). Thus, the 

w s 

variation of W and W is as is given by figure (3) . 
w s 

P(4>) must be continuous over the airfoil, so P(4>.;) “ and P(0} = PC2ir). 

JL 

“is ^ 

Substituting our values for V and a into equation (27) yields 

PC-ti) = In 2 1 cos ( I - aC<J))3 I -In (43) 

Since, at the trailing edge, P(0) = P(27r3, 

ln2|cosa^| -In jy^W^CD)W^C0)^J ° ln2|cosctjJ -InlVj Jf^(27r) 

Tit 




C413 
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Ai: all other segment boundaries, which are generally outside o£ the cusp 
region, Vi^ = = 1,P(^^") =P (4)^'*']) , and therefore 

ln2|cos[^- o.p|-ln[V^W^»iD] = ln2rcosc|i i-ln[Vj^^,lW^Win (45) 

On the lower surface, is replaced by in equation (45). 

If, in any segment, 4> - 2a^, P(4>) becomes infinite. This is most lihely 

to occur in the segments to either side of the stagnation point, 4 't 

where indicates the stagnation point. In order to prevent this, we require 

(46) 

and ' 

TT +2a, > >7r +2 a,- (47) 

L L L ^ 

X£ all the and a. *s along with Uj Uj 4>_> are given 

in the problem specification, this leaves only and left to be determined. 

However, with equations (28)- (30) there are three conditions that must be met. 
Therefore, we need to free one of the quantities listed above. The quantity best 
suited for this is 4 >t » location of the stagnation point. 

If we use equation (43) to define P((jj), equation (29) becomes 


2tt , 

/ [ln[cosC| -csCq)),|-ln V((j)) -In W^((ji)-KylnWg (4i)+ln2] coscfid^^^ir (48) 
o 

This can be expanded to 

‘^’i 

f [ln|cos(-| -a.}|-lnV^-inW^_-K^lnW^+ln2] cos(|)d4. 

.!<{>• X 

+ S ^ ^ (InJcosCl -a.)i-inV.”lnW -K InW +ln2] cos(j)d4)=7r (49) 

h-i 
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The integral cos(j) In | cos y - a. jd^ can he evaluated as 

M X 

/costjilnj cos ^ -a. [d^ = (sin^ + sin2a.) ln[cos ^ - a. | 

^ JL X 4# X 

+ 1/2 C'f’cos 2a^ - sint})) + constant 


If we denote W . as 
cx 




W^. = - f cos^ In d^ 

^i-1 


and introduce the notation 


In Ci»j) “ ln|cos - a^)| 


equation (49} can be written as 


(SO) 


C51) 


(52) 


_ la 

Kh + Kj^ W^_ + S {sin 2a^ (In Ci,i) -In (i-1, i}} 
a i=l 

1 1 

+ cos 2a. -i- -y (sin^. - sin(j)._^3+ sin(J). [ln(i,i) -InV.] 

^ X X~ X ^ X X** XX X 

dt 

-sin ^ £ln(i-l, i) -In V.]} - / ^ cos ^ JnW d^i 

X iV 

o 

2tt „ 

- / cos<{) In W C^} d(j) *ir 
^w 


(53) 


Due to equation (45} > the next to the last terra, with i=n, and the last term, 

with i=n+l (n=l, I -1}, cancel each other out. Out of these terms, only the 

last with i=l and the next to the last with i=I remain. However, since <i> =0 

& o 

and cf)i^ “ 2ir, sin = sin (|)j^ = 0, these terras also drop out. The third term 

of the summation also drops out. Therefore, equation (55) can be written as 

I, 


a 

S 

i=l 


^cl ^ Cln(i,i} -ln(i-l, i)} 


‘I*, 


w 


2 it 


i > 
) I 


r-i 


2 *‘‘^i " “^i-l^ ^“i^ ” ^ \ 'f* W (((]} d4) 

o S’ ^ 

CS4) 


= IT 
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We noiu introduce J such that 

c 


= I {sin 2«=, CClnCi.iMn Ci-l.i3) | 


1=1 


''vr 


27T 


- / COS ^ In W d({) - / cos (f> In W 

W w 




w 


I£ we define 0, we can alter the indices to get 

£L 


a 


- ir * ^ {sin lnCi,i) -sin InCi, i+1) 

i=l 




w 


+ j (}>* Ceos 2«^ - cos - / cos <jj 


2ir 

/ cos In W d<fi 
w 


<}>. 


w 


Now equation (54) can be written as 


a 


In a similar manner* if we define W . as 

S3. 

W . = / sin ({) In W C<{>) 

sr s 

'i'i-i 


and we define J„ as 
s 


Jg - X {-(l+cos 2«^) In Ci»i) + Cl***cos 


In 


i-0 




w 


+ (sin 2°:^ - sin 2«^_j,j^) } - / sin (j) In W 


27T 

“ ^ sin In W C4’) <i‘{> 
w 


0 . 


w 


-) CDS 2«.} 
“1 1 

C55) 

In d(}s 

(56) 

(57) 
CSS) 

Ci,i-J-1) 
m d<l3 

V 

C59) 
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Equation (30) can now be written as 




C60) 


From equation (45), 

In = In (i,i) -In (i,i+l) + In 

Substituting this value o£ into equation (44) yields 

lnjcos«j j -In (i,i) + In (i, i+1) + In In W^(0) 

= ln|cos«j I - In Vj - ^ In (2 tt) 


(61) 


(62) 


From equation (45), ^2 determined as a function of and so forth until 

Vj is reac±Led<, In this manner, and are eliminated from (44), yielding 

09 

-Kj, In W^(0) + ^ In W (2tt) * ^ {-In (i,i) + In (i,i+l)} = 0 

i=0 

C63) 

Defining as being the summation term of this equation, equation (63) can be 

written as -K^ In W^(0) -s- In (2Tr) ■> = 0 (64) 

We now have three equations (equations (57), (60), and (64)) for the 

three unknowns (K^^, Kjj, and tjj^ ) . Eliminating and yields 

L 


where 


7T) E 

•l " •'s °2 " 


! 

0 


(65) 

^Sl 

lnW^(2TT) + 

"si 

a 

In 

w 

CO) 

(66) 

\l 

In f^(2TT) + 

'^CI 

a 

In 

W 

w 

(0) 

C67) 

W , 

W - W 

W , 


• 


(6S) 

cl 

sl^ cl 

a a 

si 






n 


Equation (65} is a transcendental equation for * Once the value of is 
determined, and Kjj can be determined from equations (57) and (60) . 

VJith equation (45), we now have I. - 1 equations for the I values of V.. 

oL cL X 

The last equation comes from the previously unused equation (29), which guarantees 
uniform flow at infinity. 

27T ^ a ^ i IT 

/ P(i^} dtj) = ^ ^ [injcos C 

a i=.l 1 

-Kjj In - In d(^} + 27rCln 2 - In V^) = 0 

(69) 

Now that all of the and are known, P(^) can be calculated 

from equation (43), QC^i) can be calculated from equation (34), and xC^i) and 
yC4>) can be calculated from equations (32) and (33). 

Foir practical numerical calculations, the circle plane is divided into 
2N equal parts, with the positions given by 

^ (V = 0, 1, 2, ..., 2N-1). (70) 

Next, the 4’i*s (except (J). ) a.'s, K, X, p, and p are chosen. The values of V/ 

\ 1 Ci 

Wp , N , and W_ can then be calculated (equations (51) and (58)). Using 

4 ^1 

a a 

equations (66) through (68), D^, D^, and calculated, and the transcen- 

dental equation can be established. Once 4>| is determined and Kj^ can be 

It 

calculated from equations (57) and (60). PC^)) and Q(«j)} can be determined at 

each point on the circle, x(^) and y(^) are then determined, so 2N points are deter- 
mined on the airfoil. These points are equally spaced on the circle plane, but 
they are not equally spaced in the airfoil plane. 

The resulting coordinates will yield an airfoil oriented at it*s zero lift 
line. However, the angle between the zero lift line and the chord line (8) can be 
determined and subtracted from « to yield the angle of attack with respect to the 
chord line. 
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Since the values of and are determined by the closure requirements, 
there is no direct control of these values in the input specifications. 
and determine the trailing edge angle. In order to maintain some control 
over this trailing edge angle, a desired value of can be specified, 

and by iteration, varying either on the upper or lower surface (or both 
surfaces}', or K or K, the desired can be attained. 
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II. PROGRAM 

* 

The calculations required for the solution of the Eppler problem are 
carried out with the aid of an IBM-36Q/7S computer at the University of Illinois. 
The Eppler program not only determines the profile of the airfoil but also 
determines the boundary layer momentum thickness and the energy form parameter. 
However, in the present application, the boundary layer capabilities of the 
program have not been fully utilized. 

The required programs are kept in files on the PLORTS system. The 
file name of the Eppler program is EPPLER, while the file names of the re- 
quired input data are EDATA through EDATG. A sample input data deck is 
shown in Figure Al. 

The first card in an input data deck for the Eppler program is a card 
with an Alpha-numeric listing of the titles of the cards that follow. These 
titles are read in 20A4 format. It is essential that the order of the titles 
not be changed and all titles must be included on this card, even if the 
named card is not used in the program. This first card can be thought of as 
part of the program itself, as it is never changed. The remaining cards, with 
the exception of the title, are in the format (A4, 16, 14FS..2). Some of the 
data that is input through the FS.2 format is divided by a factor of 10 in 
the program, so it is important not to specify the decimal point. All the data 
should be right justified, and the program will convert the data to the correct 
multiple Ox 10, The data that is divided by 10. in the program will be identified 
in the following discussion as having a psuedo-format of F5.3. 

The manner in which the data on each card is treated is determined by the 
title, which is listed in the first 4 spaces on each card. The data is read into 



TRA1TRA2ALFAAGAMABSZ REENDEBETAPLOTTITL 
ABSZ 9200 100 

AGAM 100 100 100 100 100 100 

TITL 

U OF I HLE 1657-20-28-19 AIRFOIL 

PLOT 3470-83372128815625330711089424000 

TRA1000027 4120 4800 4290 0975 0000 4800 9200 3200 

TRA200Q027 ISO 4290 100 5000 0542 300 4000 100 8000 0499 20014500 001 

BETA -1 100 

RE 03 01633 

ENDE 
!* 


Al^ A sanjple l,nput da.ta dedc 


the program as MARKE, NUPU, and PUFF, where MARKE is the title, NUPU is an integer 
and PUFF is a 14 element array. The data is then transferred to the appropriate 
variable according to the title. 

The first title listed on the first card is the TRAl title. The TRAl 
card is the card that inputs the and oi. , The are input in terms of circle 
divisions, and the a. are input in degrees, tj)^ is determined by the program, 
so it is input as zero. The cf). and ct. are input as pairs, and up to seven pairs 
can be input on one card. If it is desired to break the circle plane into more 
than seven segments, more TRAl cards need be specified; with a maximum of four 
cards, as storage is allowed for only 28 segments. The last i{)^ must be equal 
to the number of divisions in the circle. The c{)^ must be listed in increasing 
order, including the computed value of 

Spaces 5 through 10 of the TRAl card (NUPU) are reserved for the profile 
number. If several different airfoils are developed at the same time, they 

t 

can be identified by this profile number. 

Spaces 5 through 10 of the TRA2 card are also reserved for the profile 
number, but in this case, the profile number is used only to keep track of the 
input data, as this number is not used in the program. These spaces can also 
be left blank on the TRA2 card. 

The remainder of the words on the TRA2 card define the input velocity 

function. Words 1 through 5 define the upper surface and words 6 through 10 

define the lower surface. Word 1 is <j) , given in circle divisions, and word 2 

s 

is (ji , The meaning of words 4 and 5 depends on word 3. If word 3 is 0,0, word 4 

VT 

is k and word S is p. If iTOxd 3 is 1,0, word 4 is w* and word 5 is w. If word 3 
is 2.0, v/ord 4 is p and word S is w. Words 4 and 5 are divided by 10.0 in the 
program, so the psuedo format is F5.3. 

The specification of w and w* Cwoid 3 being equal to 1,0) is recommended 
only with large values of w* , so the path of W^ is strongly curved. The process 
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* converges slowly when w' is small, and convergence is not guaranteed when p is 

negative. For less strongly curved paths, the specification of p and w is 

* . 

recommended Cword 3 equals 2.0), 

Words 6 through 10 define the lower surface in the same manner that words 1 

through 5 define the upper surface. Thus, for a symmetrical airfoil, words 6 

through 10 would repeat words 1 through 5. 

Word 11 is referred to as ITMOD, and determines the variable that is 

changed in the iteration process to set K to the specified value. If ITMOD is 

0.0, no iteration is carried out. If ITMOD is 1.0, the on the upper surface 

are altered by a factor Aa. until K attains the desired value. If Il^IOD is 2.0, 

IS 

the a. the lower surface will be altered and if ITMOD is 3.0, the a. will be 
altered on both the upper and lower surface by an equal amount. If ITMOD is 4.0, 

K is modified, if ITMOD is 5.0, ^ is modified, and if ITMOD is 6.0, K and "k are 
modified by equal amounts, ITMOD == 3,0 or 6.0 is useful for symmetrical airfoils - 
Word 12 is K , written in the psuedo-foxmat of FS.3. Word 13 is the toler- 

r 

ance acceptable in the computation, also written in the psuedo-format of FS.3. 

A suggested value for this is .001, the smallest value available in the F5.3 format. 
Word 14 is not used. 

The next card in the list is the ALFA card. This card inputs the various 
angles of attack that the pressure distribution is developed for and that are used 
in the boundary layer portion of the program. The first word after the title is 
NAL, the number of angles of attack listed, in 16 format. JiAL can be as large 
as 14. If MAL is specified as larger than 14, it is reset to 4. The next 14 £or 
.less) words are the angles of attack, in degrees, ^mritten in FS.2 Format. If 
NAL is given as a negative number, the angle of attack will be given on the 
TRAl card, where i is on the ALFA card in FS.2 format (see the sample data deck 
in Figure A1 for an example of this) . If an ALFA card is given with no angles of 
attack and NAL=0, the angles of attack of the previous profile are repeated. 


The AGAM card controls the output o£ the Eppler program. The 16 of the 
AGAM card is ignored, but 14 AGAMCi)’s are read in i-5.2 format. In general, 
the AGAM(i) ’s are either zero or not zero. If AGAMC^} is not zero, the x and y 
coordinates of the airfoil are generated. If AGAMCi) is equal to zero, only 
the transcendental equation is solved. If AGAMC2) is not equal to zero, the 
profile list will be printed, along with a velocity distribution for each angle 
of attack' on the ALFA card. If AGAMC3) is not zero, the input data and the 
solution to the transcendental equation is printed out for the initial input 
and the final iteration. If AGAM (4) is not zero, the input data and the solu- 
tion to the transcendental equation will be printed out for all iterations. 
AGAMC5) and AGAM (6) refer to the boundary layer portion of the program. If 
AGAM(5) is not zero, the program will print out a listing of the distance along 
the surface from the stagnation point, the local velocity, the energy thickness 
foxm parame-fer energy dissipation boundary layer thickness divided 

by the momentxnn thickness] , and the momentzm thickness. If AGAMCS) is equal 
to 1.0, the local Reynolds number, based on the momentum thickness and the 
local velocity is printed out instead of the momentum thiclaiess. If AGAM(G] 
is not equal to zero, the boundary layer transition point, boundary layer 
separation point, and drag (calculated by the Squire Young Method) are printed 
out. AGAM(7) through AGAM(14) are not presently used, but are reserved for 
further use. 

At the University of Illinois, most runs are made with AGAM (I) through 
AGAM (6) equal to l.Q. This results in the most complete output. An attempt 
to run with AGAM (6) equal to zero resulted in the failure of the program for 
unknown reasons. 

Card ABSZ lists the number of circle divisions, NKR, in spaces 11 through 
15. NKR must be divisible by 4, and NKR + 1 points result in the profile of 
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the airfoil. As NKR is increased, the accuracy j£ the solution increases, as 
well as the computational time required. The maximum NKR is 120, but 60 is 
usually a sufficient mrniber unless large slopes in the velocity function are 
encountered, as with a Stratford distribution. For the airfoils designed at 
the University of Illinois, an NKR of 92 was chosen. 

The ABSZ card also lists ABFA in spaces 16 through 20, which multiplies 
all values given in circle divisions. ABFA is normally equal to 1.0. It is 
necessary to change ABFA only if the number of circle divisions is changed, so 
it is not necessary to change all the input data given in circle divisions. 

If no ABSZ card is given, NKR is set to 60 and ABFA is set to 1.0. 

The RE card is used to input the Reynolds number into the program. 

Hie psuedo-format of the RE card is (A4, 6 X, 5(211, 5X, F5.3)). The first of 
the II words represents MA, which at one time was used to determine the suction 
mode. Since the capability of boundary layer suction has been removed from the 
program, this word is no longer used. The second II word is MU, the mode for 
boundary layer transition. tVhen MU is equal to 1, transition is by laminar 
separation. If , M.I is equal to 2 , transition occurs at the first decrease in 
velocity. If MU is equal to 3, transition occurs when the velocity remains 
constant throughout a step distance or decreases. If MU is 4, transition 
occurs when the natural logarithm of the local Reynolds number based on 62 
and the local velocity exceeds or equals 18.43 H ^2 - 21.74. MU = 5 is 
similar to MU - 4, except the value that In (RE) is coa^iared to is 18.43 
H 22 “ 22.10. Therefore, MU = S is a more conservative estimate for trans- 
si tion. The F5.3 word is the free stream Reynolds number, based on the chord 
length and free stream velocity. All lengths in the program are non-dimensionalized 
with respect to this chord length, and all velocities are non-dimensionalized with 
respect to this velocity. There can be up to 5 Reynolds numbers, each with its 
c\m MA and ^^U. The program will continue to read in Reynolds numbers (i:qj to 5) 


until a zero value is read as a Reynolds number. 

The ENDE card is necessary for proper termination of the program. It 
is the final data card, and indicates all data has been read in. 

The next three titles on the list are cards that have been added to 
the program at the University of Illinois, The first of these cards is the BETA 
card, which replaces the ALFA card. If a BETA card is used instead of an ALFA 
card, either a punched output is generated or data is filed into the PLORTS 
system that is used by the Stratford program. This data consists of four parts, 
written in 6F12.9 formal. The first part is DS, the increment of the surface 
distance for each x increment. There are NKR DS's generated. The other three 
parts are a velocity function (VF), and x and y coordinates of the airfoil. 

There are NKR + 1 of each of these values. The velocity function is equal to 
the local velocity divided by Cl cosot^) . The program was originally designed 
to give a punched output, but was modified to file the data directly into PLORTS. 
However, as the PLORTS system is due to be removed from the IBM-360 at the Uni- 
versity of Illinois, it will be necessary to change back to a punched output 
deck. 

The next card that has been added to the program is the PLOT card. This 
card reads data into the system that is then either punched out or filed into 
PLORTS. Nothing is done with this data by the program, as this is only a con- 
venient method of getting data into the input deck for the Stratford programs. 

The last card to be described is the TITL card. No data is on the TITL 
card, but 'this card signals that the next card is in 20A4 format, and is the 
title of the airfoil. This title will be printed in the output and inserted into 
the Stratford input deck. 

There are some restrictions on the order the cards are read in. The ABSZ 
Cif one is used), AGAM, TITL, and PLOT cards should be read into the con^mter first. 
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although not necessarily in that order. The data on these cards remains valid 
until another similar card is read into the computer. Thus, for example, if 
several profiles are to be developed with the same number of circle divisions, 
it is not necessary to repeat the ABSZ card. The next cards to be read in are 
the TRAl and TRA2 cards, in that order. Once the TRA2 card is read in, the 
profile is generated. The ALFA or BETA card is then read in, followed by the 
RE card. 'The RE card initiates the calculation of the boundary layer. If other 
profiles are desired, new TRAl and TRA2 cards can now be read in, preceded by 
new ABSZ, AGAM, PLOT, and TITL cards, as necessary. These cards can be followed 
by ALFA or BETA and RE cards if boundary layer information is desired. The ENDE 
card terminates the program after all the profiles and boundary layer calculations 
are con^lete. 

The descriptions of the output which follows assumes AGAMCIJ through AGAMC6) 
are not equal to zero. If any of these words are equal to zero, the corresponding 
portion of the output will be deleted. 

The first data listed in the output are the input data and the solution to 
the transcendental equation. This data is preceded by the title, profile number, 
iteration number, and iteration mode CO throu^ 6} . The headings of the table of 
data do not agree with the nomenclature presented in this paper. NUE represents 
the same quantity as ALPHA is a. , WS is w and F, IVHK is w' and w', DRAK is K 
and F, DRAM is p and F, ^ is and FLA is and and LAS is and Fg* 

The next data listed are the profile of the airfoil in x and y coordinates 
and the velocity distribution for each angle of attack on the ALFA or BETA card. 

At the end of this listing, the values of CM, BETA, ETA, SX, and SY are printed 
out, CM is the moment coefficient at zero lift and BETA is the angle between the 
"zero lift line and the chord line. Since all angles of attack are given in reference 
tp the zero lift line, this angle is necessary to compute the geometric angle of 


attack. ETA, SX, and SY are apparently remnants of trouble shooting the program, 
as they are not particularly useful. ETA is the nmber of points in the circle 
plane divided by the chord and tt. This term is used in non-dimensionalizing the 
chord. SY and SY are summation of the x and y coordinates of the airfoil profile. 

The last section of data is derived from the boundary layer portion of the 
program. First there are 2 tables, one for the upper surface and one for the lower 
surface. These tables list the surface coordinate, local velocity, ^2* 

If AGAMC5) is equal to 1.0, the local Reynolds number based on 5^ local 

velocity is printed in place of fig* However, nothing in the output indicates 
that this has been done, so it is important that it be noted that AGAtlCS) is 
equal to 1,0 if this data is to be used. If 3. negative number, the flew 

in the boundary layer is turbulent. 

Follovdng these two tables are listings for the upper and lower surface 
transition points, separation points, and drag coefficients. Once again, there 
is a problem of nomenclature, as the transition points are under the heading 
INS., the separation points are under the heading TRANS,, and the drag coefficient 
are under the heading SEP , . The transition and separation points are given in 
terms of surface coordinates. 


DIHEhSICN ALVCl4)«Vil4lfMjUKEKU6) 

DIMENSICK X>M3),YYy(3t 
01HE^SIC^ fie<53 tl'A(5JtHUi5 J 

DIMENSKN TSTIS;»eANT(5},.;W(4»2»141 ,SU(4,2,14) , SA(4, 2,14) 

COf'J'CW Pl{121)»XP(l2ilvYPC UUfARGt 121 J ,)«( 121) , Y{ 121) , P ( 12U , 
lPUFFa4$5ArtAMl4),nSSl2aJiiVMl21)fANI«9Q J,ALFA(90J fFKERN(30),ABSZ, 
2ABGR,l-ABGRtIB,w«,KKft,N€ HwIPfiCi JABi JST»CMfETA,ABFA,PI,8GG£NtSX, 
3DAPGtSY,IZ2fVl|TITLEtl9) 

CCKPCK /GRZX/CDK,AAC7)te(it 7) 

III I* I 

DO 4 I « 1>16 
4 HAPKEMI) ■ C. 

HAPKE « €• 

COR * .01 
PI > '3.1415<52A!;4 
BCGEK * .C114S3292S199 
10 REAC<Sil}{PAFKENU)tI«l>i{») 

1 FCRPAKIA A4) 

AES2 « 6C. 


AEFA » 1. 

9 ABGR « 36O./A0SZ 
HAeGR> .5* ABGR 
IB a .254ABSZ «-*! 
HO s 2*IB 
NKP » 2*MC 
A8S2 » KKR 


RBPKODUPffla^lTY OF tas 
fiW^AS. PAGB IS POOR 


NO » hKR +l 
CO € »*«1»IH 
AR! « *fQ ♦ a - 2*H 

H FKERMy)- ABGR*CGSGlARI*HABGR)/(SINGtARI*HABGR)*PI) 

U PE£Cf5»2)RARItE»NLFU«PUFF 

2 FGPRAT<A4bI6*14F*.2) 

DCL2 2=* 1,16 

IF (PAftKE-£C*RAPKENU) ) GJ T3 33 

12 CONTINUE 

14 WRITE*6,3J PAPKE 

3 FDPPATdX,* IKCCRPECT DATA CA^C WITH CODE. A4) 

CIC TC 11 

13 GO TC(15, 22, 333, 140,142, 25,150, 331,375,376,14, 14,14,14,14,14) ,1 

15 NtJPPC « NUPL 


J» I 

16 1*1 

17 AMU) » RUKDJPUFFm*AtiFA, ICOC.) 
XFCAM{JS.EC.C.)JST*J 

i»r*i 

ALFA(J) *PUFFm 
IFlAMUI-AfiSZ +.1* 18,21,21 

18 J*J^l 

1FU-14J 19,20,20 

19 1*1+1 

GC TC 17 

20 REAnC5,2}HARKE,NUPU,PJFF 
GC TC 16 

21 JA8*J 
ALFAlJ + U-0 
GOTO 1 I 

22 CAlt TFAFerj 
GO TC II 

25 iFfPLFFlZJ.EC.O.) GQ TC 2d 
DC 27 .1*1,5 
REPX*PUFR2=»JJ 


IF{REPX)26,28p26 
26 REUJ = l-ES’^RERX 
IFU^^INTtFUFF (2+J-IJl 
MAtJ) = IPL./100 
«U(J) = IPU/10 - 
• 27 JR ^ J 

28 CALL GPP(ALV,K‘ALfRE,Mii» JR, AAA,EBB,DDD,A2H,AZ,AKK,PE0) 

GOTC 11 

140 DO 141 1 = 1,14 

141 AGAMtI)= PLFFIIJ 
GC TC 11 

142 AESZ = POFFU) 

ABFA = PUFFI2) 

GC TC 9 

331 IPUNCH=1 
r,C TC 33 2 

333 IPUNCH=0 

332 IFiMjPU)3331 ,335,3323 

3331 NAL=-NIPL 

DO 3332 1=1, NAL 
N=IF'i:^{PLFFiin 

3332 ALVn>=ALFA(M 
GC TC 335 

3333 NAl=hUPU 
IF(MPU.GT«i4)NAL*4 
on 334 1=1, ^AL 

334 ALV(I)=PLFF ( n 

335 IFIAGAHI 2) >336,11,236 

336 N=Q 
3361 NZ=2 

IF1I2Z.GT. C) NZ = 3 
CALL ZEZni Z,M,N2T> 

WRITE(6,279J NZ.T, TITLE 

WRITE {6, 337 )NUPFC, { AL V ( Ml , M = 1 AL 1 

337 FCR^-ATtlF tyiPPCFIL ,I5,7H jaLFFA,14F7 .2 ) 

379 FCP‘/A1( A1,19A41 

IZZ = IZZ+1 
WPITE(6,330) 

338 FCRPAT152H NX Y V-D I STR I EUT ICN FOR GIVEN ALPHA) 

IFtN.NE.C) GC TC 339 

DG341 N=1,NC 

IFn2Z.GT-62) C-C TC 2361 

339 NO = N-1 
ZN= NC 

PHIH= ZN’^HABGP 
XOR = 1CC.*X(M 
VCR =. 1CC.<«V{W} 

DC 340P=1,NAL 

340 V{wj = AES{ VFCM+CCSGlPHld - AIVIHIU 
IZZ = IZZ+1 

341 WRITEt6,342 }^C,XCF^ ,YCR V CM ) ,?'=!, NALI 

342 F0PPATtI4,2FS.3,14F7.3 J 
IZZ = UZ + 1 

WRITE (6,344 )CH,OARG, ETA, jX, 

344 F0RPAT<4F CP=,F7,4,6h Bfc f a= , FS .2,9HCEG . ETA=,F5.3*4H SX=,F6*3,4H 
1Y=, F£, 3 J 

IF( lEUNCF-EC.Cl GC TO 11 

NCl=NC-l 

iiri=iiii+i 

WRITE (7,330 1 lDStn,P = l,dJll 
WPTTE17,330) t VF{r^ ) ,M= UMJ ) 


WRITE{7t33C) 

WRITE(7,330) { V(R J ,(- = l,NQi 
GC TC 11 

* 330 F0PH^TU6F12.<3n 

375 AAA=PLFFUJ/100. 

, 8BB=FljFF<2 )/100. 

D0D=FUFF(3J/100. 

AZM^FUFFt4)/lCC. 

AZ=PtFF{5)/lC0. 

AKK=PUFF(6) /lOO, 
RE0=PUFF{7}*1CCCC. 

GC TC 11 

376 READ15,37?) TITLE 
WRITEC7,378) title 

377 FCPrATUSAA) 

378 FGPN'ATl* *t1?A4) 

GO TC 11 
150 STOP 

CEB LG SLBCHK 
END 
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SUBFCUTINE ZEZ CU i IQ , IZT J 
DIMENSICN IT(4J 
DATA IT/lH^»lh 
IF( IC-2)2 ,Zt1 

1 ID = 3 
IZ = 1 
GC TC 3 

2 IZ = IZ^IC 

3 IZT = IT(IC + n 
RETLPN 

END 


SUBPCLTINE CECF(RET,UfUS»HfD2»CDiCFjH12) 

CALL FI2B(H,F12fEPST) 

IF{HJ2, 15,1 

1 CC=f 637796 1-20. 52 UUiJ*F+ 15. 70195 2) /RET 

CF=EPS7/RET 

GO TO 15 

2 H32 = -^H 

ALPFA=0 .C29^1C.93-1.95«Ai_JG10iF12))**l„7C5 

CF=ALPHA/RET=»*C.266 

T={ l.C-l.Q/H2)/SCPTiCFJ 

PI=-hl2+C2*LS/(CF+UJ 

3F(T-10.C)3,3,4 

3 Pl=-C.77380-7* (0.0 72 12- H'l 0 . C3 263-0 .C009A2'^'T M 
Gc rr 5 

4 Pl=-(iO?-16 5i-T*[2C-49i6-T*U.i^07-7*{O.C28i-O.OOQl*T) )}}/ 

1 ( 20. 5403-7* (1.1341-7^' to. 03 12-0 .0002*1) ) ) 

5 IF(PI-P1)9,6,6 

6 C3=0.Q 
lFn-10.C)8,6,7 

7 CI=1C6.2651-7*(2C.4916-T*( 1.14C7-T*IQ. 0281-0. 0001*7) ) J 
C2=ie.04Q3-T*( 1.1341-T*( 0.0 31 2-C. 0002*7) ) 

GC TC 14 

8 Cl=7.9344O+T^(C.e73OO-7*(O.395Cl-0.Cil39*T) ) 

C2=9.605 

GC TC 14 

9 P5=-3.63819+7*( C.93C59-T*t0-C9422-7*(0.0C464i-T*0.731E-04) ) ) 

IF( P!-P5)10,11,11 

10 Cl=-76.C+T+t63.7-7.C*T) 

C2=-3Q.5-(Cl+33.5)/P5 
C3-0 .0 

GC TC 14 

11 P3=-1.0722l-H*( C .0653 7 *• 7^(0. GlC46-T*t 0.000174-7*0. 6E-05) U 
IF( FT-P3)13,12,12 

12 P2 = -C.42195-T*( 0.1 73 76- r>( C.C4C22-7* (0 .001426-7*0.25 E-C4) ) J 
C3=( (6.5*P2-10.g*P3-2. 5) *t P2-P1)+(6.5^P2-2.5*P1+1.8J *{P3“P2) )/ 

1 ( ( P3-P13 *{F3-F2)»tP2-Pl> i 

C2=<6.5*P2-2.5*Pl+l.ai/iPi-P2)-C3*(Pl+P2) 

Cl=-1 .5-Fl*<C2+2.5)-C3 + Pii'Pl 
GC TC 14 

13 P4=-2.24 22 8+T*tC.453?6+T*i G.03 6 52+T*(0.CCi926-T*O.277E-C4))) 

C3=( (19 .4*P4-30.5*P5-16 .7) *( P4 -P3 1+ ( 19 .4*P4-10 .9*P3 + 10 .6 )*£ P5-P4) J 
1 (P5-P314 (F5-P4)*(P4-P3 J ) 

C2=( l9.4*P4-iC.9*P34-lC.bJ/( P3-P4)-C2*(0 3 + P4) 

Cl = -6.2-F3* tC2+10.9)-C3*P>*P3 

14 TP-C 1+C2>» PI +C3*Fr*PI*P2 

TSP= tHi2-(Hl2-l.0)*{l,0-0. 72 195/(0. 93- I- 95«AL3G 10( H12) J ) J* 

1 {jC.55216-C.B075O*H324-J.O4d55*SOFTtO.775OO*H32-l. 10667J )/ 

2 ( SORT (0.7750C*H32-1. iJt 67 )* SORT (ALPHA)*! H12* ( f-3 2- 1 3100) ) **2 1 
CC=CFfc( (TP-C. 134*7 )*(1 .lU-P I*t 1 .0 + 1 .0/H12 ) ) / 1 7SP*RET-^«*0 . 134) +F22* 

I {1.0+PI’^(1.0-1,0/H12)) J 

15 RETUPW 
E^D 
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SUBRCUTINE GPUP (H, C2, U, 1 S , CL , Z, DZ ,DOZ,V,f'Al 
CCKMCN /GRZ!</CDK,A4(7}fBd< 7) 

Z = C2*^BS(hJ 
RBT = WPE+WRE’!‘U+C2 

CALL CDCF(RET,U,LSihi02,CdTCf »H12) 

IF(«AU»lf2 

1 V = 0. 

GC TC 13 

2 GC TC{13,13U3.6,e,5,5,3),MA 

3 IF(H)4, 1,1 

4 V “ 25.*(H + 1.S8 }’»C2*LS 
GC TC 12 

5 IF(F)6,l,i 
60= eeiNAj 

PSI '= AA(MAJ 

IF lev .ME.O. )PSI = BV»AL0GUET)*PSI 
IF(PSI“1.52 )6,7,7 

7 rF(PSl-l.q9)ll,ll,G 

8 PS!=1.52 
GC TC 10 

q PSI=1.99 
113 B = C. 

11 FZ = SIGMPSI,H) 

CALL CDCF {RET,U, US »HZ,D2, CDS, CFS,H12S) 

V=(U*(CCS-< PS I + B)TCPSJ + 02 e-PSI + l-12*{B-*-PSI ) ) J/t Bi-ABStHJ-l-] 

12 IFtV)13,l,l 

13 VDIJ =V/U 
LSU =IS/L 

17 CC2 = iCF-t 2«tH12}+LSU^‘02+V0Ul*DL 
DZ = <CC - 3.*Z+USU + \fOJJ4DL 
RETURN 
END 


StaPCUTINE H2EiH,H12,EPiT i 
IF(H14,5 ,1 

1 IF(h-1.5725gi2r3t3 

2 H12 =(SQFT(H-1.515090J ) i -227 . ia220*H+724.559l6>*H-583 -601821 
1+4.02922CC 

EPST= ( (-.03172850£55*H12+»34154a5523)+H12-l-6a6094T98 )*H12 
1+2.512588652 
GQ TC 5 

3 EPST =: (H*2.221687229--^ .22625 2829) *H+1. 372390703 

H12 = (25- 71578574’^H-a9. 5321^201) + 79,87084472 

GG TC 5 

4 H 1 2= (0.C485 5-SORT ("0.77500* I- -1 .106668 ) J /I H+1 ,43 100 ) 

5 PETUPM 
END 


SUBROUTINE GRP ( ALVf NAL i .IE, MUi JF » AAA »8S8 ,OOD, A2H » AZ r AKKf REO) 

CCMWCN FlC1213,XPa2l),YP( 121) ,ARGU21),X(121),YU2U,Ptl21)ff 
lPUFFa4),AnAR(lA)*ESU203,VF(121) fAM t90J , ALFA{ 90) ,FKEPN(30 1 ,ABSZt 
2ASGP ,EAeGRt IB,MC»NKR,NC».'iUPPa, JAB,JST,CM,ETA,ABFAtPI,BOGEN»SXt 
3CARG,SY,IZZrYl,TITlEa9) 

DIMENSION ALV 1 14) , RE ( 5 ) ( £ ) , hRE ( 5 ) f H( 5 J tD2 ( 5).,CH C5 ,2 » 141 y 

1SA{5,2»14) ,5U(5,2,143 » S ARi 5 ) ,SOR ( 5 ) |CWD ( 5 } , DD( 5 1 ,RR{ 5 ) ,REI ( 5 ) 
HABLC==2 .49 
UKK=G«0 
DC 2 J=1,JR 
HiJ) « J 

RRiJ) = RE(JJ*1.E“6 
2 WP.E(J)= SQRT(PEtJ) ) 

NALA = 1 

NAtE = lAESCNAU 

IF{NAL.LT-0> NALA =t NALE 

DO 3G lA - NALA, NALE 

IF{ AGAMI5) .NE.a.) IZZ = 99 

DC 38 JL =1,2 

ND = 2=»iJL -3 

NDSD = JU - 1 

FNO = ND 

NENC = 1 + NDSO^NKft 

AP = ALVUA )/JiABGR * .5*AtSSZ * i. 

NU = INK AP+.01*FND) + NDjD 

NLNOSD NU NDSD 

DSR = DS{NUNCSD)*ABSiFL-lAT(NU)-AP} 

S=C. 

UH=C. 

IF (AGAMtSi )4,26,4 
4 IZP = IZZ + IABS{NU-NEND) -64 
IFnZP}iC,lG,6 
6 CALL ZEZnZZ,3, IZT ) 

WRlTEt6,7) IZT, title 
WPITE(6t8)NLFRCtALVCIA) 

7 F0FNA7{Al,l9A4i 

8 FOPMATUf- ,'BCUNCARY LAYER P^OF I I.E ' , 16 , 1 OX, » ALPHA=* * , F6. 2 1 

10 CALL ZEZnZZfl, IZT) 

IF{ Jl-l)16,12,l6 

12 WPIT£16,14) IZT, (PR { J) ,PlU) ,^ = i,JR) 

14 FCRMAT{ AL,13FUPPER SURF ACE , 2X , 5 ( 3X, 3HRN = ,F 6. 2, 4X,3HHU= , 1 2) I 
GCTC 20 

16 WRITEi6,18) UT, CRR ( Ji ,MUUi ,.= 1, JP) 

18 PHRMUK A1>13KCWEP SURF ACc,2X, SOX, 3HPN«,F6.2,4Xf3HMU=,I2) ) 
IFt!ZP)26,26,20 
20 CALL ZEZnzZ,l,IZT) 

WR|TE(6,Z2) IZT, (J,*j = l,JR) 

22 FGRWAT{A1,7F S ,8H U ,5ni,9H H32 ,11H C2 )) 

GC TC 26 

24 NUNDSD = NU - N'OSC 
nSR = DSiNUNCSD) 

26 ZN = NU 

CALV=CCSG{ { ZN-1 .0)*HABGrt-ALV { I A )) 

UKl=ABS< VFtNUl+CALV) 

S = S + DSR 
DC 33 J=1,JR 
OSZ = DSR 
UR = UK 

CALL GRSiH( JI,D2( J),UR, JR1,D3Z,WRE{ J) ,MU( J) ,0,O.Q,0*0,HIJ) ,C2( J) , 
IXA, XLtCCC) 

CC{JJ = C2iJ) 


IFlAGAMl5).EQ,l.)DDt J)=UKi*PR{ JJ*D2i J) 
lF(UK)3Cj28t30 
28 SARIJ)=C.C 
SUPU)=^C.O 

30 SAt?Ul= SAR{J)+XA 

33 SLR(J)= SLR!JJ+XU 
333 IFlAC-AM{5n31i35r31 

31 CALL ZEZlIZZtlilZT) 

WRIT£(6,32 )IZT,S,LK1, tH(KK) ,CC (KK) , KK=1 , JRl 

32 FOR^^AT^Ai,F7.4fF8«^,5^F10-5,F11.6J ) 

35 IFlKU“NENDJ34j36,3A 

34 UK = UKl 

IF (LKl^GT.LHK) LKK=UKi 
NU == NU *HD 
GC TC 24 

36 IF (J.EQ.l) VRITEn,376) A A A, E68 rDCD,UKK » AZH» AZ, S t AKK 
IF {J.EC-IJ ^B!TEn,377i REQ 
377 FGRMAT(F15*A ) 

DC38 J=1,JR 

CALL H12B(H( JJ ,H12F ,CTSi 

IF{F12R.GT.FA6LC)H12R = HABLC 

CWI J, JU,IAJ = {UKI*>^ (2.5<-,p>t‘H12R i ii'OZ I J)*2. 

SAtJfJUjlAJ= SARUJ 

38 SUiJ,JU,lAJ= SLR(J) 

376 FQRf^ATt8Fl0.4) 

I = 0 

IFl AGAM(61 >35,60,39 

39 IFIAGAM(5J 

40 IF {I2Z-55) A2,44,44 
42 IF(KAL‘) AArAA 

44 CALL ZE2tIZZ,3,ITZ5 
WRITE(6,7> TITLE 
WPITE(6,46)NLFRi: 

46 FCR^<AT^lF , 'SCUNCAPY LAYcR RESLLTS PROFILE', I6> 

48 CALL ZEZI1 Z2,2, ITZ) 

WRnEtfc,A )ITZ,(PE(J>,HUJ>,w = l,JP) 

49 FCR‘^AT(A1 15X,5t3FBE=,E9.^,lX,2HMU=,I2,3X)» 

IFtl.NE.C GCTC 51 

nC 60 I=RALA,KALE 
IFl IZZoGF .57 J GCTC 44 

51 CALL 2£Z(1ZZ,1,ITZ) 

WRITE(£,5C) ITZtALVlI) ,la,J=l,JBJ 

- WPITE(7,55) ALVII) 

55 FCRRAT{* SNAM ALPFA= » , F5« 2, ‘ , ^==92, 6ENO * ) 

50 FCRF‘ATLAI,6HALPHA=,F6.2,3X, 5L I 1,7H INS. ,7H TRANS. ,TH 
CALL 2E2aZ2,l,ITZ> 

WFITE16,52>!TZ, (SL{ J ,1,IJ ,SAt J ,1,1) TCW(J,l,n ,J=1, JBI 

52 FCPf'iM ( Al, 13HUPPEF SURF ACt , 2X , 5 ( 2F7 .3 ,F 7 .4) ) 

CALL 282(122,1, ITZ) 

WRITE ( 6 ,54) nZ, I SL{J ,2, II ,SA (*. ,2,1 5 ,CW 1J,2, I),J=l, JR ) 

• 54 FORMAT! Al, 13FLCWER SURF ACc , 2X , 5 ( 2F 7 . 3 ,F 7 . 4 ) ) 

DC 56 J = 1,JF 

56 CViCIJ) = C^(J,1,I) +CW(J,2,n 
CALL ZeZ (122,1, ITZ) 

VIPITFC 6,58 1 IT7, (CV^CtJ ) , J=l, JR ) 

58 rCFk'ATiAl,8t-TCTAL CC ,7X , 5( 13X , F 8 .4 ) ) 

60 CCNTINUE 
RETLPN 

Er c 


SEP. n 


SUBRCUTINE lNPtUSrU,D2rWRE iFUM) 
GD TG{1,2»3,4, 5,8,8) ,M 

1 IFaS+l.B-4) 7,8,8 

2 IFtUS-l.E-4i 7,7,8 

3 AK-21-74 
GQ TC 45 

4 AK-22,1C 

45 FUP-ALC6(iiR£’»VsPE*D2*lJ)+AK-l£-4*H 
GO TC 9 
7 FUP=22.iC 
GO TO 9 
a FUP=-22.10 
9 RETLFN 
END 


SLBFCUTINE GRS{H,C2,UK,UKi,DLtWREfHU»HAfV»Vl»HR,D2R»XAt)ajrDCC) 
CAT A EAP,EUf^,hE/.0CCC05| .C001,~i.0/ 

BI7=C. 

IFtMA )3,101,3 
lOl VI = C. 

V “ C. 

IF(hE)3,3, 1 

1 IFtAeS{UK-UKV)+ABS{UKl-UKlV)+AeS{DL-DLV)+ABS{H-HV)-.lE«6)2r3,3 

2 IF(ABStCZ*-VHPE-C2V*WREV) .GE. .lE-6) GC TO 3 
BIT = 1. 

D2E = D2E*VJP£V/WRE 

3 G2V=C2 
HV = H 
UKV- UK 
OtV OL 
WREV = WRE 
DCC - 0- 
XA = CL 

XU = DL 
XI=Dl 
VV = V 
TFIDL )A,A,5 
A FE-F 
02£= C2 
GO TC 5C 

5 LSTR =(LK1-LH)/CL 
IF(FA.LT.AJS(STP =<Vi-V)/iJL 
IFtUK 

6 IF{ (LKl-LK) /LK - 

7 D2E = .2SCC^3/{WRE*3QR7i L:»TF.M 

HE=1.61cc77 

GC TC 5C 
e IFiC2)9,S,10 

9 02E = I .66AlCfi/^PE )>»3QRT(2. ♦'DL/ aKl+LK} ) 

HE = 1.5 725 RA 
r,c 5C 
LO X = 0. 

ISTAF = C 

IF tF .17 .0. J XU = 0. 

IF(R IT)13 , 12 ,- X - I - 

Uf^ F ( L S 7 P T t r K - t B 2 y WR £th 
GQ TC A - & - 

13 IF (X -«-DL -CLV) lZil2?19 
1 ? V c = y 

UK = LK\<^-l?7F^XS 
VG= VV + VS7R+XS 
IF(H.)1A,1A, 16 
lA FAP=i.5162 

IFiHAP+eAP4h)20,20,15 

15 XA=X 

CALL I- L2E {F ,F12 ,EPSTJ 

02F = 02*ILK/LK1)**( (5.4H12H'.5J 

FE = F 

GO 7C 5C 

16 HAP = 1.515C95 

IFC HAP+EAP-FJ17, 18t la 

17 CALL UMP{LSTR,LK,C2,WRE,H7ML,FLMJ 
IF(FL'^ + ElN'}2C,16rlfi 

13 XL = > 

H=-v< 

ISTAB = C 


GO TC 12 

20 C ALL GRLPt hfCZtlJK fViREfUSrR t Dl , DZ| DDZtVG ,M A ) 

ESP = .CC2 , 

IFtKA-3 )22,22,21 

21 ESP = .OCl 

I FIX ,EC.O. JV = VG 

22 XS = XS + .54DL 
D2M = DZ + .5*002 
ZM = Z + .5+CZ 
UM=LKV+LSTR*X5 
ISG = 1 -'v . 

IFfCZM) 25,25,23 

23 IS6 = 2 
HF =' 

IFtHR- 2.130,25,25 

25 ISGV = ISG 

IF SI STAB-9) 26,45,49 

26 ISTA6 = I STAB +1 
DL = ,5*CL 

GC TC 13 

30 IFIFV - FAP+EAP)31,32,32 

31 Dl = .5*DL* S ABStH)-HAP)/tABS(h )-HHl 
GC TO 13 

32 IFIF.LT.C.IFP = -hH 
IFSHA.LT.A) VF= VV+ XS^VSTR 

CALL GRUP(hK,D2.M,UM,WRE,USTP , DL, Z« ,DZ!« ,DD2M , VM , MA) 
ISG - 3 

02E = iC2 +CC2R 
IF(D2B)25,25,33 

33 FE= (Z+DZMJ/C2E 
ISG = 4 

IF lFE-2.)34, 25,25 

34 HCTFF = FF-APSSH) 

ISG - 5 

IFl AaSlhCIFF)-.01»35,35,25 

35 HS = flEStA6SlH)-2-4ABSlHMi+FE) 

ISG = 6 

IF (FS-ESP )36,36,25 

36 IF CFE“HAF+FAci37,3e,3a 

37 DL = CL*£AestHj-FAP)/(Ad3tHi-FE) 

GC TC 13 

38 IF(P)39,39,4C 

39 HE = -HE 
GC TC 4? 

40 UE~Ll<+LSTP»CL 

call, LMP{LSTR,UE,D2E,wRE,ht ,FUFE) 
IFSFLME-EUM)42t42,41 

41 DL = rL*FUV/tFlF-FLHE) 

GC TC 13 

42 DCO = CCC + CL4VM 
424 X=X+CL 

IFSX -DL>^ 147 ,50 ,5G 
47 H = HE 
H2 = 02E 

GO TC (43,13 ,^3,45,13 ,43 J ,ISG\» 

43 IF(HS-.1*ESP)A4,13,13 

44 DL - 2.*CL 
XSTAE =TSTAB-1 
GC TC 13 

45 IFtHCTFF)44, 13, 13 


49 STP-0.0 

TF( ISTAB.NE.9) STP-FLJATA i/ I ISTAB-9 J ) 

ViRnE(6,5U 

51 FrP*^AT(lXt*’*****WAFMNG'l4^=*‘* E.L. INFCP»ATI!3N IS f'lOT VALID AFTER 
IFERE. step size IS AT ITS L3VJEST BCUNO.*) 

50 HR = HE 
D2B = D2E 
RETURN 
END 


I 


SUBPGUTINE CRflW[WCtWS»WLtL)KTC« ,FL,AG,HA) 

CCSP = CCSG(AG*FL) 

C CALL ERAW(^C>VS,ULfDRAK( JJ J),FLA{ J),ABGR,NKPJ 

C VvE^^ NKP=OtHRD CK ALS kappa ALFGEFASSTtSGNST ALS K. 
FK= CK 
IF(MA)2tir2 

1 FK= FK*U.-CCSP)/n.+CGSPi 

2 SIMP = SING(AG*FLI 
V^L = -D^*ALCG^ I.+FKi 
IF(FL,EC.O.J FK = 1* 

BETA =(CCSP-1.)/FK + CQSP 
BQ^l = BETA**2 - 1. 

WUBEQ= SCFT lABSCBQMli J 

U- a« + eETA)=tSINP/U,4-G0SP) 

IFS BCMl )3,5,4 

3 WF = ALCGiAESiCWLBEO+UJ/lWUBEC-U))) 

GG TG 6 

4 WF = 2.* ATAMU/KLBEQi 
GC TC 6 

5 WF = C- 

6 WC -{WUBEO’^WF - SINP - BcrA*FL*AG^‘l .?4532g£-2)*D« 

= tCrSP-1 .) + tl.-(l./FK + 1. J-^ALaGU.+FKJli'DM 
RETLFM 
END 





SUBPCUTINE TRAPRQ 

Dlf'EKSICN PLReS(i3) ,FLSi2) ,FL^ { 2 J ,DRAK ( 2 J » ORAMC 2 ) , AC ( 4, 3 ) , D( 3} , 
IWSI (2) ,WCH 2), PINT (3) , A t4J jFK(2),P{ 3) 

CD^*^C^J Pl(12n,XP{12iJ,YP( 12 1) fARGt 12 1) 9 XI 12 1 ) » Y 1 121 ) tP (121 ) f 
lPUFF(14)9AGAM14J,DS(120J,VF(121i,ANI(90)f ALFA(SOJ,FKERN(30J ,ABSZ, 
ZASGRtFAECtR rIG♦^'C,^Kft,^lC^NJPfiG, JAB»JST9CN»ETA,ABFA,PI ,BOGEN,SXt 
3DARGtSVf!ZZ tVl, TITLE (19 J 
III = 99 
DO 23 I = 1,13 

23 PUPES^IJ=PL^D(PLFF^I), 1000.1 
1 = 1 

J=1 

24 FLS(J)= PLPESd J=*A£FA 

CALL i:P.A(s(hC,feS,tiL ,. 6 ,-l.,FLS( J),ABGP,U 

CALL' CRA^^WCIlJ)T^^S^^J),rfl.I ,-.6r”l- »FLS( JJ ,A8GR,1) 

WCK J J= V»CI ( J )+UC 
USI(J)= l^SI(JJ+(iS 
WLI = WLI+WL 

FLA(J1= FUPESn-Hl *■ ABFA 
IF(FLA( JU23,25,26 

25 CRAK(J)= 0 
DRAM(J)= 1. 

GCTC 34 

26 WI = CCSG( ABGR*FLA{J>) 

IF (PUBES £142 )-l .127,30,29 

27 OPAKm= .L=tPUPES(I+3J 

28 CPAf^(J)- .l’«PtPFSa4-41 
GCTC 34 

29 C9AK(J) = (( .I4PLDES(i+43 -lO./PUPES ( I + 3J J-1. ) *ll-+WI J/(lo-V>I ) 
DRflK(J|)= RUNC(DRAK{J),10 jO.J 

CRAViJl = -!=*PUFES(i-»-3> 

GC TC 34 

30 AA = .C5»(1.-V>I )+PUKESU + 31 
WIL^ = ALOP {.l*PUFESil44n 
FPIT = ,5 

N'lT = C 

31 FW = -WTLN/iiLCG(AA/FMIT +1.J 
>^IT=*^IT4 1 

IF( ABS(FP-FMT 1-1. E-61 33,32,32 

32 FMIT = FN' 

GC TC 31 

33 CRA^ ( JJ = BL■^C{F^,lCOO.i 

DPA5 = .C5*PLPES( 1*3) 4( Wi«-1 . J/FW 
DPAKU) = i^L^C(DRAS,loa0.i 

34 1= 1+5 
J= J + 1 

IF (J-3) 24 ,38 ,38 
39 WER = 0 

U SI t 21= -WSl 12) 

Gil) = BLi ^ (VSH2 J+-i^S H 1) ) 

0(21 =-m *{k^ci(2)+wciai i 

C(3) = KCK (2)-hCI(2l»kiSI(i ) 

ITMCU = FLF.ES(ll) 

7T>-R = IT?/cd 

SHKS = . 1>*PLFES( 12) 

HKST = .!< PL^^ES ( 13) 

35 DC 3e J=l,3 

36 ACd , J)= 0. 

ALI\< = C. 

SINAI = C. 

Cr’SAI = 1. 


FNI = 0, 

J = 1 

37 CSAIP = CCSG{2.*ALFAC J) i 
SNAIP s SING(2.+ALFA{J) > 

IF(J“JST“1140,39,40 

39 ACiZrU= SINAI 
AC{2,2i= ~1.~C0SAI 
AC (2,3)=^ -1 o 
ACU,1)= -SNiSIP 
ACO*2J= 1.4CSAIP 
AC(3,3i“ 1. 

AC(4, 1J= CCSAI-CSAIP 
AC(4t2J= SINAI-SNAIP 
AC(4,3J= C* 

ALIS- == ALIV 
ALISF = At,FA(JS 
GGTO 41 

40 FII = CSLGCHAeGR*FM-9C.,ALIVJ 
FIIP= CSLC-(HAeGF+FMi-90.fALFAt J) J 
PB-FM*hAeGR=teCGEN 

AC{ 1 fU = "FI IF^SNAIP*fiI^$iNAH<COSAI-CSAIPl*PB -i-ACdf U 

AC( lf2)=-FI I*a.+CGSAiJ+Fi IPi«C l.+CSAIPJ 4{SINAI-SNAIP)*PB+ACtlt2J 

AC(1|3)= Flip - FII + ACtIt3# 

41 IF( J'vlAE“l)42»43t43 

42 ALIV =ALFA.(J) 

SINAI=SNAiP 
CnSAI=C5A!P 
FNI = AMU) 

J=J + l 

GC TC 37 

43 DC 4? J=l,2 
iF(f=LA{JJ) 47*47,49 

49 CALL CPAW{WC,VkS,i*L ,CRAKt JJ *D3AN(J},FLAtJ )*A6 GRtOI 
ACU,1)= WCi AC(1,1) 

IFU-2)45r44,45 

44 V.S = -VtS 
WL =•' WL 

49 ACa,2) = WS *4C(1,2) 

ACUrlJ =-WL 4AC(1,3) 

47 CGNTIK^uE 
DC 92 J=l,4 
A(J) = C. 

DC 92 1=1,3 

£2 A( JJ = A( J)- + nt I )*AC U,I i 
LGTPA 
£3 r = c 

FV = 9.B9 

PHIS'H = .5 +ULIS4ALI SP) 

60 CSLI = CSLG IFHISF ,ALiSi 
CSLIf= CSLGCFHISHjALISP) 

FP=(A(U+A{4)*90 .♦EOGEN ) A t2 ) *CS LI +A 13 I^CSL I P ) +AC 4) I SH^BGGEN 
lf{ I"20}£iTe6,66 

61 IF4 AEEC FPl-AES (FV ) + .££-9 i 6 2 , < 3, 63 
63 ! = 2C 

PHISh = FHISF - PCIF 
GC T'' 6C 

62 PCIF = -FP / lA 12}/tPriIi.H-ALfS } + A ( 3) /( PHISH-ALISPH 
I = ! + i 

65 FV=FF 

PFISF = PHISF + PDIF 
GC 6C 



I 


66 ANIiJSTJ = (PhISh+^0. J/HAdGR 

70 DC 71 1=1,3 

71 FINK I J-AC( 1,TJ +AC(2,I3^'Ci>LI + AC{3,I )*CSLIP4-AC{4,I)*BCIGEN + {PHISH+ 
F90. ) 

69 FHtl J = {FIM{l)*WLT-FINT(33*ViC U2J J/Ct2J 
HK(2J=( FIM (l)*^iLHFINT(3J*WCI( 1) )/0I2) 

HKS = HK(l)+FKi2l 

IF{/iE5(FKS -SHKS) -HKST ) 7^ , 72 , 7 2 

72 IF( ITPCC J73 ,75,73 

73 IF(AGiMtA) 376,98, 76 

irr^cn =c 

75 IFlAG/lPO) J 76,3CC,76 

76 NZ = 2 

IF{ I2Z+JAB.CE.59) NZ-3 
CALL ZEZU 2Z,NZ,NZT) 

WPITE(6,7)NZT, TITLE 
89 WRITE(6,77)NUPPC,MER,ITi'lK 
7 F0Rh'^T(Zl, 19/14) 

77 FCRM^TUF ,20HRESLLTS FGil PPCFILE , 15, 5X, lOHITERATIQN , I3,5X,5FPC0E 
1 ,12) 

IZZ = IZZ+1 
fc,P!TE?6,7ej 

78 FCRPATC* KLE ALPHA nS )*HK DRAK DRAM HK FLA LAS 

1* ) 

Jh= 1 
JN= 1 

79 IZZ = IZZ+1 

IF { J^-l )fC, ei, 8C 

80 IF ( JN“J/ie)82,Hl,g3 

81 XI = ',5»U.+ C^SG (FLA< JHJ-^ABGR J ) 

WhK = {l. + CF^K(JH)’^{i„-XA)/Xl)^*(-OFAMlJH) ) 

KSTR= CPAR< JF)^CFAK IJH )/Xi 

WPITEt6,82) /INK JM),ALFA IJW ) , W S7R,WHK , DR AK { J H ) , DRAM! JH) , FK t JH ) ,FLA C 
IJH ) ,FLS{ JH5 

62 FCR‘'A1{F7.?,F6.2TF7.3,Fi:.i,2F7 .3, F9 .6, F5 .1, F4.1 j 
JF= JF+1 

IFt JF-J^E)65 t £6,85 
83 WF ITEt6,62 ) AM( JN) ,ALFA lJi< ) 
f'5 JN=JN + 1 
GCTC79 

66 IFt IINDOICC ,?,G0,100 
SB !F(MFt; ) ICC,99, 100 
9Q I F( AGAP{ 33 ) 76, ICC, 76 

100 IF { TTf^nn-4J lOl, 120,120 

101 IF(f^FP } lC3rlC2,lG3 
10? DAL = .1 

G'^ TT 104 

103 DAL '= { SFKS“HKS)4CAL/(HKS-hKSV ) 

CAL = FL'NCt CAL, 100.) 

IF(CAL)1C4,7^,104 
' 104 J = 1 
HL = C. 

105 IF{hUlC7,lC6,lC7 
'106 IF{ ITMr.D-2) IC6, 1C9, ICa 
IC7 IF tT'^Hno-i ) ica, in,iod 

108 ALFAIJ) = ALFA( J)+CAL+tl.-HL) 

109 IF(J-JST)lll,llO,lU 
lie FL=2, 

111 J = J+I 

IF (J-J A 6-1) 105,112,105 

112 HKSV=FKS 


MEK *MER4-1 
GCTC 35 

120 IF(NER)122»121,122 

121 DDK=.l 
GCTC 123 

122 DDK = (ShKS“>HKS)*CDK/(hKS”riKSV) 

D0K= PUNC(DCK?1000.) 

IF(DCK)123t 7'«f,123 

123 If { I7MCD-5) 124,125,124 

124 DRAK(1)= CF!AK(1)+DCK 

125 IF {ITMQC-4U26, 112,126 

126 DRAK(2)= OP JK <2 ) -*-CCK 
GC TC 112 

300 IF(^GAMII) )301, 11,301 

301 AK = - .5*ICCSG«PHISH-ALFAUST + 1) )/SlNGlPKISH-ALFA( JST + IJ ) 

1 -CCSGIPHISH-ALFAi JoT J )/ SIMGC PFISh-ALFA( JST ) ) } 

AKP =AK+18G./S. 8656044 
PHIM = C. 

NU = 1 
1 = 1 
ANU s=C, 

JH = C 
VI= C. 

302 JH=JH+l 

FFl = CCSG( AS3R*FLA( Jh) J 
FF2 = DRAKUH)/{1«*FF11 
FGl = CCSGt ABGR*FLS( JH) i 
FG3 = ,6/(FGl-l.) 

304 \l- VI - CGLGI PHIV-90. , ALFAtD) 

GC TC 310 
306 APGl - ANU 

IFI ANU.GT..5* ABSZ3APGA= ABSZ “ ANU 
CSP = CCSGUPGM+A6GR1 
F = n . 

IF (AH GK.LT. FLA tJH) )F=CRA>lJMi'^LCGC iCSP“FFU ’*«FF2 + 1 - J 
G = 0 . 

IF(AFGK .LT.FLS (JH) }Gs-hK tJ H ) ♦A LCG ( 1 ( CSP-FGl ) *FG3 S * *2 ) 

P(*NU}= F + G+CSLGlANU^^HABGrl-gC. , ALFA{ I) J +VI 

PKM ) = FfMU )-AK*ABS(SING( ( aNU*FABGR ~ 9C.1 - PHISHJ) 

NU = MU + 1 
A^U= AM+ 1, 

310 !F(ANU-AM( IM306, 306,3 12 
312 IFIAM'- ABSZ )314 ,320,220 
314 PH!H = AMt I )^HABGP 

VI = VI + CSLG( FHP-90 .,ALFAtn ) 

I = I+l 

IF( I- 1-JST) 3C4, 302,30^ 

320 PS=C. 

B2=C. 

DC 324 !=lT^KP 
Fs=ps+p m 
BI = 2*( I“1 } 

324 B2 = B2 + S I NG ( Bl’t A3GR J 1 1 ) 

VI = 2.»EXP IPS/AeSZ) 

SXI = .CCCCCCCC 
SY=C. 

00328 N=l,NC 

Q = 0. 

IJC326 H = l,r8 

WN = N + 1 + MO - 2'i'M 


j.'-* 



IF .GT.NKP ) NN = MN - iNKR 
IFIWN.LT*!) = Rh- + NKR 
326 0 = 0+ FKERN(N)*{Pl(HN)-PilFWn 
ANU= N-1 

Z° = ANL^HABGR - 9G. 

ZL = COSG(ZP - PHISHJ 
ZU = ABS( (1 --ZLI/Cl.+ZU) 

IFIZL .NE.O. )2L=A10GUU 

ARG^^) = 0 - AKP*SINGCZP-PHI3F)+ZL + 2P 

VFIN; = V1*EXF(-F(N) 1 

WV = CnSGIZP l/VF{N J 

XP(N)= V(V*S ING (ARG4NU 

yP{N)=-k*V*CCSGCABG<Nl ) 

sxr= SXI+ XP(NJ 

328 SY =' SY + YFfN) 

SX = SXI 

XPK = SX/lAfcSZ - 1.) 

329 X(1)=C. 

YPK = £Y/(AESZ -!•) 

SXI = .CCCCCCOO 
YU ) = C. 

XPUJ= XPK 
YP{1)= VPK 
XP{ NC)=XPK 
YP^^C J = YFK 
HCV = 0* 

C033CN = 2rNC 

SXI= XP(N-11 -XPK - XPK + XP(N) + SXI 
X(N) = SXI 

V(N} = YlN-1) + (YPIN-U - YFK + YPIK) - YPK) 

RQ = X(^)»X(N )+Y{MJ*Y(N) 

IF(Pf:,GT.R9VJL=N 

330 RQV = RC 

DC 327 I = 1,3 
lEPFL = L-2+I 

32 7 R ( I )=SQP7{X (lEPPU *XUEP^u)+Y( lEPFL lEPPL) J 
333 TA'J = (o(3)-P{ UJ/{4.*(R(2J+Ri 2)-P( 1)-R(3J) ) 

XF.AE = X(U +TAU*(X{L + i)-XI L-IJ +2 .*T AU* ( X (L+-1 ) +X ( L~1 J -Xt U-X ( U ) ) 

VNAR = Y(U+TAU*lYa + l)-YlL"'l)+2.‘‘TAL*(Y{L+l)+Y{L-lJ-YtU-Y{Ul 3 

SO = XNAS-XPAS + Y^AS»Y^Ai 

AT=XNA5/Sn 

«= YNAS/SQ 

STREF - l./SCRTtSC) 

ETA = ABSZ + S7P.EF/PI 
CV - ,9=<ETAiSTREF*B2 

CARG = 1S.C9E59 ^ ( 3 ,*YM AS/ XMA 5 - t YNAS/ XNAS 1 +*3 ) 

SX = £'rReF+SXt200. 

SY=STREF*SY^2C0. 
on 331 ^=2,^C 
XR = X I N J 

X(M= l.-B«Y (N )“AT*XR 
Y(N)= B*XR -AT+=Y(N) 

ARG(N) = ARG(M} - CARG 

WO = ( XP;^l) tXPfN-l 3-XPK-XPKJ +*2 + £ YP ( N) +VP (N^U -YPK-YPK 1*«'2 

331 DS(^-1) = STPEF*SCRT£ WQJ +£ l-+*£666667*( tXPCN3*YP(N-l} 
L-XP{^-1)*VP^^ J )/WG}**2) 

X(l) = 1. 

332 ARC- { 1 ) = ARC f 1) - CARG 
U PETUFM 

E^D 



SUBfCUTlNE CIA(X,V,hPrDJ 

DIHb^SIC^ U(3)rV{31,Wl3) ,Zl3),/i(3)tBC3i fX(121J fY(121) 

D = Q* 

DO N=2,NP 
NR ^ NP-1 

1 IF! XtNP l-X{N))3,2,2 

2 NR - NR - 1 
GO TC X . 

3 DNN=Y{N)-Y(NR)-{YtNK+l)-YiNR) M (X(N }-XI NR) * /tXtNR+lJ-X(NR) ) 
IF(CNN-C}5,A,4 

^ D ^ CNN 

5 DG 6 I = 1 * 3 
NO = N+l-I 
m - NR-i+I 
U(I)=X(NC) 

VXI}=Y(NC) 
i^(I )^X(NU) 

6 ZiIJ=Y(NU 

7 CALL CtPtlJ,V,A) 

GALL CIFCVJfZfE) 

IF{ABStA(3)-B(3J) .LT..COJJ,) GO TO 11 
XST = (Bt2)-Al2))4.5/CAtii-Ei3J) 

YS = A(2)+2.*AnJ*XST 
IF{ AES(YSi-.CCCUS t8,8 
R TA - l./SCRTCl«+YS=*YiJ 
DC 10 I = 1*3 
XQ = TA*tL( U + vm*YSJ 
vm = TA-^( vn i“YS=*u{Xi } 

U(I) = xc 

XO = D-fYS^ZilJ ) 

zin = zn )-YS=*w( I) ) 

10 wen = xc 

GC TG 7 

9 0^ H+iM2)-Et2> + ( A( 3)»e{ 3n’*XSTJ*XST 

11 ^ETUFN 
END 


FUNCTIC^ SINGIA) 

SING =OSIN( A*. 0174532925 199DC0 J 

RETUI5N 

END 


FUNC7IGN COSGtA} 

CCSG =DCCS{ A+.D1745329251S9DCG ) 

RETLBN 

END 


FUNCTIC^ PUND(i»ej 

RUNG = { AINT(A^0+SIGN( .SjAi n/B 

RETUPN 

END 


FUNC7ICN CSLG<ArB) 

CSLG = AUOGiiBSJSINGlA-ajJ J 

RETIPN 

END 


